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Abstract-  Fourier  transform  methods  are  the  standard  way 
for  determining  time-domain  pulse  structure  and  arrival 
time  from  a  set  of  continuous  wave  (discrete  frequency) 
underwater  acoustic  model  calculations.  This  technique 
requires  a  large  number  of  computer  model  runs  at  closely 
spaced  frequencies,  often  making  it  computationally 
expensive.  It  has  the  advantages  of  including  the  correct 
attenuation  at  each  frequency  component,  and  of  correctly 
treating  continuity  requirements  at  the  water/sediment 
interface.  Direct  time-domain  computer  models  are  not  as 
accurate  for  ocean  bottoms  with  strong  attenuation  over  a 
large  bandwidth  of  frequencies.  In  this  work  the 
frequency-domain/Fourier  approach  is  optimized  for 
maximum  efficiency  at  a  given  level  of  acceptable 
imprecision.  Techniques  are  presented  to  improve  the 
efficiency  of  the  individual  frequency  component 
calculations,  and  to  avoid  running  many  of  the  frequencies. 
Efficiencies  at  individual  frequencies  are  gained  through 
intelligent  selection  of  grid  parameters  in  the  ocean 
acoustic  model  (a  parabolic  equation  model).  Further 
improvements  are  achieved  through  intelligent  zero 
padding  schemes,  and  by  interpolating  envelope  functions 
at  the  receiver  location  in  order  to  estimate  (and  hence 
avoid  running)  up  to  90%  of  the  calculations  required  by 
the  Nyquist  sampling  theorem.  The  effects  of  the  various 
approximations  are  shown  in  the  examples. 

I.  Introduction 

Fourier  synthesis  of  time  domain  (TD)  results  from 
continuous  wave  (CW)  model  calculations  is  the  standard[l] 
way  of  accurately  modeling  pulse  propagation  in  underwater 
acoustics.  Direct  TD  modeling  is  possible[2],  however  an 
accurate  treatment  of  frequency-dependent  bottom  attenuation 
over  large  bandwidths  remains  a  challenge.  For  accurate  TD 
results  the  procedure  involves  making  many  CW  runs  over  the 
relevant  frequency  band,  each  correctly  treating  the  frequency- 
dependent  bottom  attenuation,  followed  by  Fourier  synthesis. 

The  practical  limitations  of  synthesizing  TD  results  from 
many  CW  model  results  are  largely  imposed  by  the  total  time 
window  of  the  result  (7)),  which  is  governed  by  the  frequency 
spacing  of  the  CW  runs  (Af)\  Ti  =  1/A f 


While  it  is  desirable  to  include  the  entire  time  for 
propagating  a  signal  pulse  from  the  source  to  receiver  point,  in 
practice  this  may  necessitate  a  very  small  value  for  Af,  which 
will  consequently  require  a  large  transform  size  with  a 
correspondingly  large  number  of  CW  model  runs. 

A  simple  optimization  to  somewhat  alleviate  this  frequency 
sampling  problem  is  to  zero-pad  the  frequency  band  outside  of 
the  source’s  significant  bandwidth,  thus  eliminating  many  of 
the  otherwise  required  model  runs.  While  certainly  better  than 
running  all  frequencies  and  subsequently  multiplying  results 
outside  the  source’s  bandwidth  by  zero,  even  this  optimization 
seldom  improves  runtimes  by  more  than  a  factor  of  two. 

Another  commonly  used  optimization  is  to  choose  a  time 
window  long  enough  to  contain  the  pulse  at  the  receiver 
location,  but  not  necessarily  long  enough  to  contain  the  entire 
propagation  time.  The  total  propagation  time  can  later  be 
estimated  from  range  and  a  reference  sound  speed  in  the  water 
column,  and  the  total  propagation  time  obtained  by  adding  on 
subsequent  time  window  values  until  the  estimated  time  is 
approximately  reached.  There  are  two  potential  problems 
associated  with  this  method:  (1)  the  estimated  propagation 
time  may  not  be  accurate  in  complicated  underwater 
environments,  thus  making  the  stacking  of  time  window 
lengths  imprecise;  and,  (2)  the  time  window  chosen  may  not 
be  long  enough  to  contain  all  multipath  arrivals  of  the  pulse, 
thus,  leading  to  wrap  around  in  the  final  TD  result. 

The  optimization  method  developed  here  is  to  make  a 
limited  number  of  model  runs  at  regular  spacings  across  the 
bandwidth  of  interest,  and  to  then  interpolate  the  pressure,  as 
amplitudes  and  phases,  to  a  sufficiently  high  sampling  rate. 
This  results  in  significantly  reduced  runtimes  for  the 
underlying  CW  model,  with  little  degradation  in  the  final  TD 
result.  In  subsequent  sections  of  this  work,  this  method  is 
detailed  in  the  context  of  a  developmental  example,  and 
directions  for  further  research  are  described. 

II.  Method  and  Developmental  Example 

The  method  described  here  and  the  physical  motivation 
behind  it  are  best  illustrated  with  the  aid  of  a  elementary 
example.  For  simplicity  and  manageable  computational  times, 
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a  Gaussian  pulse  centered  at  100  Hz.  with  a  width  of  5  Hz. 
was  used  as  the  source.  The  test  environment,  including 
source  and  receiver  locations,  is  shown  in  Fig.  1.  The  long 
propagation  range  (141  km)  and  seamount  near  90  km  were 
chosen  to  give  a  variety  of  very  difficult  challenges  to  the 
methods  tried.  The  parabolic  equation  underwater  acoustic 
propagation  model,  RAM [3],  was  used  to  simulate  the 
acoustic  propagation.  A  transmission  loss  (TL)  field  plot  at  the 
center  frequency  was  generated  and  is  shown  in  Fig.  2.  A  TL 
level  of  roughly  -95  dB  for  the  center  frequency  at  the  receiver 
point  shows  that  there  will  be  significant  acoustic  energy 
present  in  the  water  column,  even  at  this  relatively  distant 
range. 


0  50  100  150 

Range  (km) 


Figure  1 .  Test  environment  for  developmental  example. 
Source  position  is  shown  on  origin  axis,  receiver 
POSITION  IS  SHOWN  AT  141  KM  RANGE. 
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Figure  2.  Transmission  loss  field  plot  at  the  center 
FREQUENCY  OF  100  HZ. 

A  reference  TD  solution  was  generated  using  conventional 
Fourier  synthesis.  The  RAM  model  was  used  to  obtain  the 
CW  acoustic  propagation;  512  individual  frequency  runs  were 


made  in  a  band  50  to  150  Hz.  Zero  padding,  as  discussed 
above,  was  applied  to  frequencies  where  the  source  weighting 
function  was  more  than  100  dB  down  from  the  center 
frequency  strength;  this  corresponded  to  100±24  Hz.  Out  of 
the  original  512  frequencies,  only  246  actually  required  model 
runs,  requiring  300  seconds  on  a  desktop  computer1.  Because 
of  the  sparse  frequency  sampling  rate,  the  propagation  time  to 
the  141  km  receiver  range  was  wrapped  an  estimated  18  times. 
The  resulting  output  pulse  is  shown  in  Fig.  3. 
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Figure  3.  Fourier  synthesized  output  pulse  from  test 

SOURCE  AND  ENVIRONMENT. 

Examination  of  the  complex  pressure  components  of  the 
pulse  as  a  function  of  frequency  (Fig.  4)  shows  a  sinusoidal- 
like  regularity.  A  simple  way  to  reliably  interpolate  between 
the  calculated  frequencies  is  not  evident.  For  clarity,  this  and 
subsequent  figures  will  only  show  the  band  100±5  Hz, 
however  the  relevant  features  are  representative  across  the 
entire  bandwidth. 
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Figure  4.  Complex  pressure  components  near  the 

CENTER  FREQUENCY.  SAMPLE  POINTS  ARE  SHOWN. 

1  2  GHz  Macintosh  G5  computer,  using  the  OSX  10.4 
operating  system  and  Absoft  Fortran  7.0  compiler. 


Transforming  the  complex  pressures  into  their  phase  and 
magnitude  components  greatly  simplifies  the  structure,  as  is 
shown  in  Fig.  5.  The  magnitude  varies  slowly  and  regularly  in 
a  narrow  band  at  roughly  ±20%  of  its  average  value,  while  the 
phase  appears  to  be  an  approximately  linear  function  wrapped 
in  the  interval  from  -n  to  ±7C.  This  and  several  other  trial 
cases  have  shown  that  the  phase  varies  regularly  in  the  manner 
seen  here  with  a  fairly  consistent  slope  that  may  be  positive  or 
negative,  depending  on  the  underwater  environment  and 
receiver  placement.  Given  the  goal  of  eliminating  the 
necessity  of  running  many  of  the  intermediate  frequencies,  the 
regularity  shown  in  these  figures  makes  it  clear  that  at  least 
some  success  may  be  had  using  interpolation  if  a  way  to 
reliably  unwrap  the  phase  can  be  found. 
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Figure  5.  Phase  (a)  and  magnitude  (b)  of  the  complex 

PRESSURE  NEAR  THE  CENTER  FREQUENCY.  SAMPLED  POINTS 
ARE  SHOWN. 

If  it  can  be  guaranteed  that  the  phase  wrap  count  can  be 
accurately  maintained  between  sample  frequencies,  the 
frequency  domain  may  be  sampled  more  coarsely  than 
indicated  in  the  Introduction.  Note  that  this  will  not  violate 
the  Nyquist  sampling  requirement^]  for  transforming 
between  the  frequency  and  time  domains,  because  it  only 
requires  that  the  complex  pressures  be  sampled  at  or  above  a 
certain  level.  There  is  no  requirement  on  the  model  (or  method) 
by  which  the  pressure  values  are  generated. 

One  very  simple  algorithm  for  unwrapping  the  phase 
reliably  for  large  frequency  strides  is  to  thoroughly  sample 
near  the  center  frequency,  and  to  obtain  an  average  frequency 
slope  in  this  region.  Trial  unwraps  in  more  sparsely  sampled 
frequency  ranges  can  then  be  made  by  linearly  extrapolating 


along  this  slope  out  to  the  next  phase  point,  and  counting  the 
number  of  unwraps  that  will  bring  the  extrapolated  phase 
nearest  to  the  target  phase.  By  keeping  track  of  the  total 
unwrap  at  each  calculated  frequency,  the  phase  relationship 
shown  in  Fig.  5a  can  be  reduced  to  a  nearly  linear  function. 
Subsequent  interpolation  to  the  desired  sample  rate  can  then 
be  made,  followed  by  transformation  into  the  final  TD  result. 

The  results  of  this  procedure  applied  to  the  example  case 
are  shown  in  Figs.  6  and  7.  The  frequency  stride  for  the 
interpolated  result  is  1.8  Hz.,  or  10  times  the  original  sampling 
rate.  The  section  ±2  Hz  around  the  100  Hz  center  frequency 
has  been  fully  sampled  to  obtain  the  average  frequency  slope 
for  the  above-described  unwrapping  procedure.  The  overall 
speed  gain  from  this  optimization  is  a  factor  of  7.3. 


Figure  6.  Interpolated  and  original  phases  in  the  (a) 
85-95,  (b)  95-105  (corresponding  to  Fig.  5),  and  (c)  105— 
115  HZ  FREQUENCY  BANDS.  MlS -WRAPS  ARE  SEEN  NEAR  94 
AND  104  Hz. 


Figure  7.  Interpolated  and  original  magnitude  of  the 

PRESSURE,  ACROSS  THE  ENTIRE  FREQUENCY  BAND. 

A  number  of  mis-wraps  in  the  interpolated  phase  caused  by 
the  increased  frequency  stride  are  evident  in  Fig.  6,  and 
several  large  deviations  from  the  original  pressure  amplitude 
are  seen  in  Fig.  7.  However,  when  transformed  into  the  time 
domain,  the  final  result  is  in  good  agreement  with  the  original 
solution,  as  shown  in  Fig.  8. 


Figure  8.  Time  domain  results  from  the  phase  and 

MAGNITUDE  INTERPOLATION  OF  THE  CASE  DESCRIBED  IN  THE 
TEXT,  OVERLAYED  ON  THE  ORIGINAL  RESULTS  FROM  FlG.  3. 

The  RAM  model  was  used  to  generate  the  examples 
presented  in  this  paper.  However,  the  technique  presented, 
which  converts  multiple  CW  results  into  broad  band  pulse 
results,  is  applicable  to  any  CW  propagation  model.  A  ten-fold 
decrease  in  overall  computational  time  is  typical.  It  is  possible 
to  achieve  significantly  greater  reductions  in  overall 
computational  time  if  the  technique  is  mated  to  a  particular 
propagation  model  and  both  technique  and  model  are 
optimized  for  the  ocean  environment  under  consideration.  As 
an  example,  consider  that  the  arrival  times  and  general  pulse 


structure  are  desired  for  an  acoustic  pulse  traveling  in  a 
shallow-water  waveguide  environment.  The  center  frequency 
of  the  pulse  is  1  kHz  and  the  acoustic  bandwidth  is  1  kHz, 
extending  from  500  Hz  to  1,500  Hz.  It  is  possible  to  choose 
different  propagation  parameters  for  the  RAM  model  that  are 
critically  sensitive  for  each  CW  frequency  component  used  in 
the  Fourier  synthesis;  similarly,  constraints  can  be  lessen  on 
less  sensitive  parameters.  When  combined  with  the  technique 
discussed  above,  the  result  shown  in  Fig.  9  is  obtained. 


Figure  9.  Time  domain  results  for  a  1  kHz  bandwidth 

WITH  CENTER  FREQUENCY  AT  1  KHZ.  THE  FULLY  SYNTHESIZED 

PULSE  IS  SHOWN  IN  BLACK  AND  REQUIRED  300  MINUTES  TO 
GENERATE.  THE  PHASE -MAGNITUDE  INTERPOLATION  RESULT 
SHOWN  IN  RED  TOOK  40  SECONDS. 

The  black  curve  is  the  result  obtained  by  running  the  RAM 
model  for  each  frequency  over  the  entire  1  kHz  bandwidth  and 
performing  a  Fourier  synthesis.  The  RAM  model  parameters 
were  standard  selections  that  ensure  high-fidelity  propagation 
predictions  for  each  particular  acoustic  frequency.  The  red 
curve  is  the  result  obtained  by  applying  the  technique 
discussed  in  this  paper  together  with  RAM  model  parameters 
that  were  chosen  to  retain  only  the  essential  physics  needed  to 
generate  a  broad  band  result  from  sparsely  selected  CW 
results.  The  fully  synthesized  pulse  (black  curve)  required  300 
minutes  on  a  desktop  computer.  The  sparse  synthesized  pulse 
(red  curve)  was  obtained  in  only  40  seconds  on  the  same 
computer.  While  not  a  perfect  replica,  the  red  curve  retains  the 
correct  first  arrival  time,  very  nearly  the  same  second  arrival 
time,  and  preserves  the  overall  structure  of  the  pulses.  There  is 
a  450  times  reduction  in  run  time.  A  future  paper  will  discuss 
in  detail  how  the  phase-magnitude  interpolation  method  can 
be  combined  with  the  RAM  propagation  model  to  produce 
rapid  broad  band  simulations. 

III.  Future  Development 

In  the  work  presented,  a  basic  observation  of  the  behavior 
of  acoustic  phase  and  amplitude  with  respect  to  frequency  in 


Summary 


an  underwater  waveguide  environment  has  led  to  the 
development  of  a  technique  that  can  rapidly  produce  a  TD 
result  from  a  sparse  number  of  CW  calculations.  The  initial 
results  are  very  encouraging  and  give  insight  as  to  where 
improvements  could  be  made. 

Further  research  is  needed  in  the  phase  unwrapping  part  of 
the  algorithm.  Using  variances  in  the  slope  of  the  phase  in  the 
sampled  band,  or  sampling  the  average  slope  at  several  places 
in  the  source  band,  could  produce  a  more  reliable  unwrapping. 
A  more  robust  approach  is  to  apply  a  series  expansion  to  the 
total  pressure  field.  Accurate  phase  estimates  may  be  possible 
by  using  information  obtained  from  the  derivative  of  the 
magnitude  function. 

Another  improvement  under  consideration  is  a  more 
sophisticated  interpolation  scheme.  In  the  work  presented,  a 
linear  scheme  was  used.  Higher-order  interpolation  algorithms 
could  improve  the  magnitude  interpolation. 

Finally,  because  a  TD  pulse  is  being  reconstituted  from 
heavily  decimated  information,  this  method  may  be  useful  as  a 
sound  compression  algorithm.  The  requirement  of  starting 
with  complex  pressures,  and  possibly  the  inability  to  encode  a 
continuous  signal  as  opposed  to  the  isolated  pulse  in  this 
example,  may  prove  to  be  practical  limitations  on  the 
method’s  usefulness  in  sound  compression  applications.  This 
interesting  ancillary  application  warrants  further  investigation. 


IV. 

A  rapid  method  for  computing  and  transforming  underwater 
acoustic  CW  model  results  into  TD  signals  has  been  presented. 
The  method  is  based  on  the  observed  behavior  of  pressure 
phase  and  magnitude  in  an  underwater  waveguide 
environment.  An  example  of  pulse  propagation  using  this 
method  was  presented  with  a  net  speedup  of  7.3.  When  the 
method  was  combined  with  environmentally  optimized 
propagation  model  parameters,  and  allowances  for  acceptable 
error  were  considered,  a  speedup  of  450  was  realized. 
Directions  for  possible  future  work  were  discussed. 
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